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ABSTRACT 

We examine the effects that the modelling of a Boxy/Peanut (B/P) bulge will 
have on the estimates of the stellar gravitational potential, forces, orbital structure 
and bar strength of barred galaxies. We present a method for obtaining the potential 
of disc galaxies from surface density images, assuming a vertical density distribution 
(height function), which is let to vary with position, thus enabling it to represent the 
geometry of a B/P. We construct a B/P height function after the results from a high- 
resolution, Wbody+SPH simulation of an isolated galaxy and compare the resulting 
dynamical model to those obtained with the commonly used, position-independent 
“flat” height functions. We show that methods that do not allow for a B/P can induce 
errors in the forces in the bar region of up to 40% and demonstrate that this has a 
significant impact on the orbital structure of the model, which in turn determines its 
kinematics and morphology. Furthermore, we show that the bar strength is reduced 
in the presence of a B/P. We conclude that neglecting the vertical extent of a B/P 
can introduce considerable errors in the dynamical modelling. We also examine the 
errors introduced in the model due to uncertainties in the parameters of the B/P and 
show that even for generous but realistic values of the uncertainties, the error will be 
noticeably less than that of not modelling a B/P bulge at all. 
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1 INTRODUCTION 

Many edge-on disc galaxies can be seen to contain boxy-, 
peanut- or X-shaped isophotes, which are usually grouped 
together into the category of Boxy/Peanut (hereafter B/P) 
bulges. As a result of a number of theoretical studies (see 


ject), including orbital structure and stability analysis ([Bin- 

ney| 

1981[ |Pfenniger||1984| |1985| Patsis et al.||2002| [Skokos 

et a 

l.||2002a|b[ Martinez-Valpuesta et al.||20061 [Harsoula & 

Kalapotharakos||2009| |Contopoulos & Harsoula||2013 ), and 

numerical simulations ( 

Combes & Sanders] 1 1981 

Combes 

et al. |1990| Raha et al. 

|1991| |Mihos et al.| 

|1995 

Athanas- 

soula & Misiriotis|2002| 

Athanassoula|2003[ | 

12005 

IBureau & 

Athanassoula |2005| [Martinez- Valpuesta et al. |2006|), these 


structures are now known to be due to vertical instabilities 
in the bar, which cause it to ‘puff up’, giving rise to boxy or 
peanut-like shapes. These studies also show that once a bar 
forms, a B/P bulge will form soon after. 

Observational studies have further confirmed the link 


between B/P bulges and bars (see Kormendy & Kennicutt 
2004 and Kormendy 2015 for reviews on the subject), by 


showing that the fraction of edge-on disc galaxies with B/P 
bulges is comparable to the fraction of disc galaxies con¬ 
taining bars ( Liitticke, Dettmar Pohlen|20QQ ). Kinematic 


studies of edge-on barred galaxies and B/P bulges also con- 


firm the connection between the two structures (Athanas- 

soula & Bureau 

1999] Bureau & Athanassoula| 1999] [Chung 

& Bureau] |2004 

[Bureau & Athanassoula[[2005[ and refer- 


ences therein). Therefore, barred galaxies at present and 
past epochs will contain B/P bulges, and in fact, one is also 


believed to be present in our own Galaxy (Weiland et al. 



Howard et al.[[2009| [Shen et al.[ 

2010] [McWilliam & 

Zocca 

112010] Ness et al. 2012, 2013a,b; 

Vasquez et al.|2013[ 

Wegg & Gerhard[[2013| [Gardner et al. 2014[ [Nataf et al.[ 

2013] 2014 2015). 

) thirds of disc galaxies in 
strengths ([Eskridge et al. 

Bars are found in about twc 
the local universe, with variable 

2000| [Menendez-Delmestre et al. 

2007] [Barazza et al.[[2008[ 

Aguerri et al.|2009||Gadotti|2009 

), and are known to be the 
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main drivers of the secular phase of galaxy evolution. The 
torque they induce into the disc causes outward angular mo¬ 
mentum transfer, which in turn will cause a redistribution 
of matter in the disc. They are thus responsible for driving 
gas to the centre of their host galaxy ( Athanassoula| 1992b ), 


forming discy pseudo-bulges ( Kormendy fc Kennicutt|2004 























































































2 Fragkoudi et al 


Athanassoula|2005 ), redistributing stars in the galactic disc 
(Sellwood & Binney 20Q2| Roskar et al. 2008 Minchev & 


Famaey 2010), and possibly creating a fuel reservoir for AGN 


activity ([Shlosman et al.[1990 iCoelho fc Gadotti|2011 Em- 

[sellem et al.|2014| but see also |Lee et al.|2012[ For a review 

see |Gombes|[200l| . However, even though the effect of bars 
on all these processes has been thoroughly examined, a study 
of the effects of B/P bulges on all these processes has not, 
until present, been carried out. 

As a first step towards understanding the effect B/P 
bulge geometry may have on the aforementioned processes, 
in this paper we focus on their infiuence on dynamical mod¬ 
els of their host galaxies. These models are obtained directly 
from images of the galaxies’ surface brightness, by first as¬ 
suming a vertical density distribution, or height function, 
and subsequently deriving the potential of the galaxy. They 
have been used extensively in the literature, with one of 
their most important implementations being in simulations 
that study the response of gas in a fixed potential. These re¬ 
sponse simulations are used to study the dark matter content 


of galaxies and to test the maximum disc hypothesis (Kranz 


et al.|20dT] I Weiner et al.|2001||Slyz et al.|2003[|Kranz et al. 


2003[ [Perez et al. 2004), the bar pattern speed (Lindblad, 
Lindblad Athanassoul^|1996[ [Kalapotharakos, Patsis 


Grosb0l|2OlOb| ) as well as the kinematical and morphological 


properties of gas in galaxies ( Lin et al.|2011 2013). Dynam¬ 
ical models have also been used in studies determining the 
bar strength ( Buta fc Block|2001| Laurikainen fc Salo|2002| 
and the orbital structure of galaxies (Quillen et al.[|19^ 


Patsis et al.|1997| [Kalapotharakos et al.|2010a[ [Patsis et al.] 

2010| . Furthermore, they have been used to study gravita¬ 
tional torques in barred and spiral galaxies in order to estab¬ 
lish the amount of gas infiow and by extension determine the 
importance of secular evolution ( Zaritsky fc Lo||1986 Haan 


et al.|200^ Foyle et al.|2010 ). In all of these aforementioned 
studies, the geometry of the B/P bulge is not taken into ac¬ 
count when constructing the height function, and instead a 
position independent, ‘fiat’ height function is assumed. This 
is partly due to the lack of an analytical model for a B/P 
bulge, as well as to the inherent difficulty of detecting these 
bulges in face-on or intermediate inclination galaxies, which 
are the galaxies generally used in these studies. 

Various methods however have been proposed over the 
past few years, which allow either for the detection of B/P 
bulges, or at least for an educated guess at their existence. 
By viewing a large number of iV-body+SPH simulations, 
and covering a wide range of viewing angles, |Athanassoula| 
et al. (2014) have shown that B/P bulges manifest them¬ 


selves in face-on projections as the so called ‘barlens’ (Lau¬ 


rikainen et al. 2011), which renders their detection fairly 


easy. Strong observational arguments for this have been pre¬ 


sented in Laurikainen et al. (2014). Another method pro¬ 


posed by Debattista et al. (2005), uses signatures in the stel¬ 


lar kinematics of face-on or almost face-on galaxies and was 


implemented by Mendez-Abreu et al. (2008) who confirmed 
the existence of a B/P bulge in NGG 98. Furthermore, it is 
possible to detect signatures of B/P bulges by examining the 


morphological features of inclined galaxies ( [Athanassoula fc 
Beaton||2006 Erwin V Debattista||2013 ) 


We therefore believe that a study of the effect of B/P 
bulges on models of their host galaxies, and by extension of 
the necessity of including the geometry of B/P bulges in the 


height function of these models, is called for. To this aim, 
we first introduce in Section a straightforward and reli¬ 
able method for calculating the potential, forces and deriva¬ 
tives of forces of a general density distribution p(r,0,z). We 
present some tests which demonstrate that the method can 
give highly accurate results, while also allowing the fiexi- 
bility to choose an arbitrary height function, without being 
restricted to one which is constant with position. 

We then used our code on an image of an V-body+SPH 
simulated galaxy, which is presented in Section [231 fhus ob¬ 
taining a realistic potential for a barred galaxy. In order to 
create this model, we assign a thickness and a height func¬ 
tion to the galaxy. These height functions are introduced 
in Section and include two ‘fiat’ height functions, a func¬ 
tion which describes peanut bulges (from which we construct 
our fiducial B/P model), and another which describes boxy 
bulges. 

The main results are presented in Section where we 
examine the effect B/P bulges have on the potential and 
forces (4.1), on the periodic orbits (4.2) and on the bar 
strength (4.3). We find that B/P bulges indeed have a signif¬ 
icant effect on the results and therefore conclude that they 
should be included when modelling their host galaxy. 

In Section!^ we explore the errors which will be induced 
in the results by using a B/P model which is not exactly 
the ‘correct’ one. This is necessary since it is not trivial to 
observationally obtain the exact parameters of B/P bulges, 
and this can introduce errors in the model. We show that 
for a range of uncertainties, the errors induced in the results 
are less than those induced by not modelling the B/P at 
all. We also introduce a new method for calculating the bar 
strength, which takes into account the variation of the 

non-axisymmetric forcings along the whole extent of the bar. 

In Section we give a summary and list the main con¬ 
clusions of our work. 


2 METHOD TESTS 

To create a dynamical model of a galaxy we first need its 
density distribution. The two-dimensional surface density 
can be obtained from surface brightness images of a face-on 
disc galaxy by assuming a M/L ratio. It is important that 
these images are taken in a wavelength range which min¬ 
imises the effect of dust extinction and traces the old stellar 
population (for example, the Spitzer 3.6/xm band). In this 
work we use an image of a face-on simulated galaxy, and thus 
we do not need to account for dust extinction, nor assign a 
M/L ratio, as our two-dimensional image gives the surface 
density directly. Once we have the surface density we also 
need to assign a height function, and together these give us 
the three-dimensional density distribution, from which we 
can calculate the potential of the galaxy due to the stel¬ 
lar component. Our method for calculating the potential in¬ 
volves a straightforward three-dimensional integration over 
the density distribution and we refer to it throughout the 
paper as the 3DF method. We calculate the potential in 
Gartesian coordinates by 
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X (kpc) 

Figure 1. Comparison of orbits in the analytic and 3DF poten¬ 
tials of our model galaxy. The thin black line gives the outline of 
the bar. The orbits calculated in the 3DF potential are given in 
solid thick black lines, and the orbits calculated in the analytic 
potential are given in dashed red lines. We plot some xi orbits 
(along the bar), some X 2 orbits (perpendicular to the bar), and 
an almost circular orbit outside the bar region. We see that the 
orbits in the two potentials almost completely overlap, so that 
the red and black lines are practically indistinguishable. 


^{x,y,z) = 

/ oo POO no 

-oo J —oo J —a 


-G 


fMd^£^=dx'dy'dz\ ( 1 ) 


- Xj) 


+ e2 


where G is Newton’s gravitational constant, p is the density 
and e is the softening length which is necessary to elimi¬ 
nate the noise at the expense of a small bias ( Merritt [1996 


Athanassonla et ST]|2000 ). We can differentiate the expres- 


sion in Eq. analytically with respect to x, y and to ob¬ 
tain expressions for the the force and its derivatives. We thus 
rely heavily on an adequate integration algorithm, specifi¬ 
cally one which can deal with singularities. To tackle this 
we use CQUAD, a doubly adaptive integration algorithm 
( Galassi et al.|[2003 ), which requires more function evalua¬ 
tions than other integration routines, but is more successful 
in dealing with difficult integrands. It computes the integral 
within the desired relative error limits (or precision), which 
the user can set. Since we mainly work in the z=0 plane, 
we focus in what follows on the non-zero quantities in the 
plane: the potential $ and the two non-zero components of 
the force Fx and Fy. As mentioned, the above, as well as 
what follows, concerns the potential and forces of the stellar 
component of the galaxy. 


2.1 Tests on the method: Relative errors 

In order to test the accuracy of our method, we create a 
model of a barred galaxy containing a disc, a bar and a clas¬ 


sical bulge, using density distributions which have analytic 
solutions for the potential and forces. We then calculate the 
potential and forces for this model using the 3DF method, 
and compare the results against the analytic solutions. The 
general results of these tests are very positive, which demon¬ 
strates the ability of our code to deal with difficulties such as 
cuspy and/or non-axisymmetic density configurations. For 
more details on the model of the galaxy we refer the reader 
to Appendix [A| 

To calculate the relative errors of the potential and the 
derivatives of the force, we use the relation 


Error = 0.5 


\Ri-R2\ 

\Rl\F\R2\' 


( 2 ) 


where Ri and R 2 are respectively the analytic and 3DF so¬ 
lutions. 

To calculate the relative errors for the force, we use the 
relation 

\Fii — Fi2\ . . 

Error = ' A, (3) 

where i and j can be either the x or the y component of 
the force, and the subscripts I and 2 stand for the analytic 
and 3DE solutions respectively. We therefore normalise the 
error of each component of the force by the total force at 
each point. This is done because our main interest in the 
forces is for the calculation of orbits and because on the 
symmetry axes of the x- and y- components of the force (in 
the static frame of reference), the analytic estimates of Fx 
and Fy will be equal to zero. 

We stress that the precision with which the code cal¬ 
culates the results is an input parameter to the code. The 
accuracy can be as high as the user wants it to be (within 
the limits of machine precision), at the expense, of course, of 
computation time. We require a three dimensional integra¬ 
tion and, due to the propagation of error at each integration, 
the relative precision we ask of the CQUAD algorithm for 
each integration has to be larger than that which we wish 
to achieve. Practically this means that if we ask for a rel¬ 
ative precision of I0~^, we will obtain a relative precision 
of approximately 10“^. This is sufficient for our purposes 
as the error is less than 1% for all variables. The softening 
is set to 10 pc throughout the paper. Eor this precision and 
softening, the maximum error of the potential is 0.3%, of Fx 
0.6% and of Fy 0.7%. 


2.2 Tests on the method: Orbits 


Even though from the results of the relative errors we see 
that the 3DE method gives highly accurate results, we would 
like to confirm that the noise in the force field does not 
prevent orbits from running smoothly, as they would in an 
analytic potential. To do this, we calculate a number of pe¬ 
riodic orbits in the analytic potential and in the potential 
derived using the 3DE method for the galaxy in the frame 
co-rotating with the bar, in the model described in Appendix 
\A\ The grid used for the orbits and throughout the paper is 
200 X 200, and the orbits are calculated using a Kick-Drift- 


Kick leapfrog algorithm (Hockney & Eastwood 1 1988 Quinn 


et al.|I9^ Springel||2005 ) 

In Eig. [2 we plot a few of these orbits. In this figure 
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and all throughout the paper the bar major axis is along the 
x-axis. In the figure we show some xi orbits, which are elon¬ 
gated along the bar, some X 2 orbits which are perpendicular 
to the bar, as well as some nearly circular orbits outside the 
bar region. We see that the orbits calculated in the 3DF po¬ 
tential are a very good approximation of those calculated in 
the analytic potential, as the two practically coincide. Thus 
the error which is introduced in the potential from our 3DF 
method and the adopted value of precision is sufficient for 
our purposes. 


2.3 The image 


We use the 3DF method on the density distribution derived 
from a face-on image of a simulated isolated galaxy and dif¬ 
ferent height functions (which are described in Section . 
The initial conditions of the simulation from which the im¬ 
age was constructed, include a live spherical dark matter 
halo, an exponential stellar disc and 75% gas (for more in¬ 
formation on the simulation the reader is referred to run 
gtrllG in Athanassoula, Machado Rodionov| ( |2013 )). The 
snapshot we use is taken well into the secular evolution phase 
of the galaxy, specifically at 8 Gyr after the start of the sim¬ 
ulation, by which point a strong bar and B/P bulge have 
formed. The image we use is constructed from the ‘stars’ 
component of the snapshot and has a morphology reminis¬ 
cent of that of many strongly barred galaxies, such as IC 
4290 (see Fig. |^. 

In order to decrease the noise of the image we require 
a snapshot with a large number of particles. We create a 
snapshot with 40 times the particles of the original snapshot, 
following the procedure described in Athanassoula (2005). 


To further reduce the noise in the image we apply some 
smoothing, by Fourier decomposing and recomposing it. 
The Fourier components are calculated as follows: 


dnir) = — /* S(r, ^)cos(n^) d^, (4) 

^ J — TV 


bnir) = — /* ^{r,0)sm{n0) dO, 
^ J — TV 


(5) 


0 is the azimuthal angle, r the radius and S gives the sur¬ 
face density. We then reduce the high frequency noise by 
recomposing the image as 


= Y + XI {(^m{r)cos{m0) + bm{r)sm{m0 )), (6) 


using only a limited number of even Fourier components (in 
our case nF=26). We show in Figs. 3(a) and 3(b)| the surface 
density of the original image and of the Fourier recomposed 
image, respectively, both in arbitrary units, and in Fig. |3(c)| 
we show the residual image of the two. As the images of 
surface density are in arbitrary units, the density, as well 
as the potential and its derivatives will also be in arbitrary 
units in what follows. 


3 HEIGHT FUNCTIONS USED 

In order to obtain the three-dimensional density of a galaxy 
disc from a two-dimensional image we need to assume a 
height function, which defines how the density drops off as 
a function of ^ from the equatorial plane z=0. The height 
function and the scaleheight (zq) will of course affect the 
results, and we therefore need to use the height function 
which best approximates that of the galaxy we are trying to 
model. 

The height function can be either constant or can 
change with position. In the case where it is constant with 
respect to position we assume, for simplicity, that the den¬ 
sity distribution can be written as 

p{x, y, z) = 'Six, y)F{z), (7) 

where p is the three-dimensional density distribution, E is 
the two-dimensional surface density, and F is the height 
function. In the more general case where the height function 
depends on position in the galaxy, as would be for example 
the height function describing a B/P bulge, the scaleheight 
changes as a function of position. In this case, the density 
distribution would be given by, 

p{x, y, z) = S(a:, y)F(x, y, z), (8) 


where an and bn are the even and odd Fourier components. 




(a) IC 4290 


(b) gtrllG 


Figure 2. Visual comparisson between IC 4290 and gtrllG. The 
two galaxies have striking morphological similarities and are clas¬ 
sified as having the same bar strength (for more details see Section 

iFSl. 


where the normalisation of the height function is 

/ oo 

F{x,y,z)dz = 1. (9) 

-OO 


It is worth noting that the mass of the model is always 
the same; the height function simply determines the volume 
density of the galaxy, by setting the thickness of the disc. 


3.1 Flat height functions 

Up to now in the literature, position-independent or ‘flat’ 
height functions have been used when modelling barred disc 
galaxies. We therefore also use two flat height functions in 
this paper, to check the discrepancy which will be created in 
the model by a) using a flat height function and a B/P height 
function, and b) using two different flat height functions. We 
adopt two commonly used functions, the isothermal-sheet 
model ( van der Kruit fc Sear le|l 981 ): 
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Figure 3. Left: Original image showing the surface density of the stellar component of the gtrllG simulation. Middle: Model from the 
Fourier recomposition, using up to 71^=26 even Fourier components. Right: Residual image after subtracting the model from the original 
image. 


F[z) — ^;^sech^ 
2^0 



and the sech-law model ( van der Kmit||1988 ): 


( 10 ) 


F{z) = —sech ( —) , (11) 

TTZo \zq J 

where l/(2zo) and l/(7rzo) are the respective normalisation 
factors. 


3.2 Peanut height function 

To obtain a height function for the peanut, we examined 
the particle distribution along different cuts in x and y from 
the simulation introduced in the previous section. We found 
that the sum of two two-dimensional gaussians for the scale- 
height can provide a reasonable approximation to the B/P 
shape. As can be seen in Fig. and as commented below, 
this choice may fail at certain points, but provides an overall 
fair representation of the structure. 

The resulting B/P height function is a non-separable 
function of position and is given by: 

F{x, y, z) = — 1 — rsech^ ( f / ) . (12) 

22 o(a:,y) \zo{x,y)) 

The scaleheight 2(0 (x,y) varies like the sum of two two- 
dimensional gaussians: 


a(i,») =^exp (- )) + 

(13) 

where A is the maximum scaleheight of the peanut above the 
disc scaleheight, The variance of the gaussians is given 

by cr^, (xo, yo) is the position of the maximum of the first 
gaussian and (xi, yi) the position of the maximum of the 
second gaussian. We fit these two two-dimensional gaussians 
to values of the scaleheight obtained from the simulation 


along y = 0 and x = 3 (which is where the maximum of the 
scaleheight occurs). In the remainder of the paper, we refer 
to this as our fiducial peanut (or B/P) model. 

To obtain the scaleheights, we take cuts along the x- 
and y- axes and fit the vertical particle distribution with 
a sech^ function. We thereby determine the variation of zo 
from bin to bin along the cut. The results can be seen in 
Fig.[^ In the side-on view (panel (b)) we see that the scale- 
height along y = 0 behaves approximately like the combi¬ 
nation of two gaussians, except in the central region where 
the scaleheight drops below that of the outer disc. For a cut 
along X = 3, where the peanut is maximum (end-on view, 
panel (c)), the behaviour of zo is still well approximated by 
a gaussian, although our fit slightly under predicts the value 
of the scaleheight. 

Along some cuts at x values intermediate between the 
centre and the peanut maximum, the gaussian approxima¬ 
tion fails to represent the behaviour of the scaleheight with 
y. In fact the behaviour of the scaleheight in the presence 
of a B/P bulge is quite complex, and cannot be grasped en¬ 
tirely by a simple analytic function. However, as it turns out, 
the fitted function shown in Fig. [^underestimates the value 
of Zo at these points. This directly translates into an under¬ 
estimation of the effect of the peanut in those regions. In 
summary, our fiducial model for the peanut height function 
shown in Fig. |4]will result into a conservative estimate of the 
effect of the real peanut present in the image we adopt as 
our starting point. Given the scope of this paper, which is to 
demonstrate the generic effect of a peanut bulge on a galaxy 
model, we find this approximation more than satisfactory. 


3.3 Boxy height function 


The B/P bulge might at times have rather boxy isophotes. 
This could be due to projection effects, whereby the peanut 
is projected at such an angle that the isophotes appear boxy 
( Athanassoula Misiriotis||2002 ). However, boxy isophotes 
might be present even when the bar is seen side-on, i.e. they 


might be the real shape of the B/P bulge (see Patsis et al. 


(2002) for a discussion based on orbits). This tends to be 
the case for galaxies with weak bars, where instead of a 
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0.0033 0.023 0.1 0.42 1.7 

(a) gtrll6: Side-on 




(b) Edge-on: Side-on 


(c) Edge-on: End-on 


Figure 4. Left: side-on image of the surface stellar density of the simulated galaxy gtrllG. Middle: The scaleheight of the simulation 
(red crosses) is plotted along the a^-axis (for y=0, i.e. the side-on projection). The solid black line shows the fit of 2;o(ic,0) to the data, 
which gives the scaleheight of the fiducial peanut height function. Right: Plot of the scaleheight of the simulation (red crosses) along 
x=3 which is where the maximum of the peanut occurs (end-on projection). The solid black line shows the values of zo{x,y) along x=3. 


strong x-shape or peanut forming, boxy isophotes are seen 
( Athanassoula||2006 ). 

To model a boxy bulge we use a height function which 
drops off as sech^ with height from the z=0 plane, 


F{x,y,z) = 


2zo{x,y) 


sech 


zo{x,y) 


(14) 


where the scaleheight is a top-hat function. 


( Zq ^ \x\ ^ Xmax ^ \y\ ^ Vruax 

otherwise. 


(15) 


4.1 B/P effect on potential and forces 

We calculate the potential and the forces for the density 
distribution given by the image described in Section |2.3 
and the different height functions described in Section [3 
The results of this subsection are shown in Fig. In the 
top row we plot the scaleheight along the bar major axis for 
the two setups we are comparing in the plots. We show two- 
dimensional plots of the relative difference for the potential 
in the second row, of the x component of the force in the 
third row, and of the y component of the force in the fourth 
row. The green line represents an ellipse fitted to the outer 
isophote of the bar. Each column gives one of the following 
comparisons (from left to right): 


and where Zq'^^^ gives the scaleheight of the disc and 
gives the scaleheight, or strength, of the boxy bulge. This is 
quite a simplified model of the boxy bulge, with only two 
free parameters, its strength and length (which is set by 
L=2xmax)- The thickness of the box, i.e. ymax, is set by the 
width of the bar. 

We create the fiducial boxy height function such that 
it best approximates the fiducial peanut height function, so 
such that we can examine whether the former can be used as 
an approximation for the latter, as there is one less parame¬ 
ter to model. The fiducial boxy height function therefore has 
a height equal to the height of the fiducial peanut model and 
its length is such that the boxy bulge finishes approximately 
where the peanut scaleheight is in between its maximum and 
minimum (see top right panel of Fig. [^. 


4 BOXY/PEANUT OR NO BOXY/PEANUT? 


We wish to investigate whether accounting for the geometry 
of the B/P in the height function will significantly change 
the model of its host galaxy, and therefore whether we should 
include it in the modelling when it is present. Thus, in the 
next three subsections we investigate the effect B/Ps will 
have on the potential and forces (Section 4.1), on the peri¬ 
odic orbits (Section |4.2| and on the bar strength (Section 


4.3) of the model. 


(i) Two flat height functions: sech and sechf 

We compare the models obtained by implementing two 
flat height functions, sech and sech^, with equivalent scale- 
heights, in order to demonstrate that different flat height 
functions do not significantly affect the results. In the very 
centre, the difference for the potential is only around 1% and 
for the forces it is 5%, while for the rest of the grid the dif¬ 
ference between the two setups is well below 1% in all cases. 
If we decrease the value of the scaleheight the two height 
functions produce even more similar results. This happens 
because as the disc tends to become infinitesimally thin, the 
shape of the height function becomes less important. Equiv¬ 
alently, if we increase the value of the scaleheight, and hence 
the thickness of the disc, then the difference in the results 
obtained with the two height functions increases. 

We see therefore that the scaleheight, and not the vertical 
profile of the height function, is primarily responsible for 
creating differences in the models. 

(ii) Flat and peanut height functions 

We compare a flat height function and the height function 
of our fiducial peanut bulge, i.e. a peanut with parameters 
fitted to our simulated galaxy. The differences that arise 
from using these two setups are significant for the potential 
and forces, as can be seen in the second column of Fig. 
This is especially true near and around the region of the 
maximum height of the peanut, and in general in and around 
the bar. The force can be different in the two cases by up to 
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Figure 5. Errors from not taking into account the proper B/P geometry: The top row gives the scaleheight (in red and black 
lines) along the bar major axis for the setups we are comparing in the plots. The second, third and fourth rows give the relative difference 
between the two setups being compared for the potential, Fx and Fy, respectively (see the colourbar for values of the relative difference). 
The dark green line represents the ellipse fitted to the outer isophote of the bar. First Column: Difference between the sech and sech^ 
setup. Second Column: Difference between the fiducial peanut height function and a sech^ height function. Third Column: Difference 
between a boxy height function and a sech^ height function. Fourth Column: Difference between a boxy height function and our fiducial 
peanut height function. We see that not including a peanut or a boxy bulge where there is one will induce large errors in the potential 
and forces and also that a boxy height function is not a good approximation for a peanut height function. For details see the text (Section 
|4TI. 
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40%, which is not surprising since around the maximum of 
the peanut the scaleheight is more than three times the value 
of the scaleheight of the disc. This can be seen in the top row 
of the figure, which demonstrates how the scaleheight varies 
along X for the two height functions. The larger scaleheight 
reduces the forces in the plane of the disc, due to a reduction 
of the density in the plane. 

Therefore, we see that by not taking into account the 
geometry of the B/P bulge, we induce significant errors in 
the model, i.e. in the potential and its derivatives. 


(iii) Flat and boxy height functions 

In the third column of Fig. we compare a flat height 
function and the height function of the fiducial boxy bulge. 
We see that a boxy height function will also induce large 
differences compared to the flat height function, and in 
fact in the central regions our fiducial boxy height function 
has an even larger effect than the peanut. Boxy bulges are 
usually associated to weaker bars in simulations and are 
therefore typically less strong than peanut bulges (although 
at early times boxy bulges can be as strong as peanut 
bulges - see Figs. 2 and 3 in |Athanassoula fc MisirioB^ 
(2002)). In observations, boxy bulges can appear as strong 
as peanut bulges (see for example Chung & Bureau (2004)) 
although it is hard to distinguish whether these are truly 
boxy bulges, or simply peanut bulges seen at an angle. It is 
therefore reasonable to assume that a substantial amount of 
boxy bulges will be somewhat less strong than our fiducial 
boxy height function. However, for the sake of simplicity 
and to be able to compare our fiducial boxy height function 
with our fiducial peanut height function, we give the 
former the same strength as the latter. Therefore our 
fiducial boxy height function can be thought of as an upper 
limit for the effect of a boxy bulge on the model of its galaxy. 


(iv) Peanut and boxy height functions 
We compare our fiducial peanut height function to our 
fiducial boxy height function in order to see to what extent 
the boxy height function approximates the peanut height 
function as it has one less free parameter than the peanut 
height function. These results can be seen in the fourth 
column of Fig. For the potential and forces the match 
between the boxy height function and the peanut height 
function is quite poor, especially in the central region 
where the scaleheights of the two height functions are very 
different. Therefore the boxy height function is not a good 
approximation to a peanut bulge. 


4.2 B/P effect on periodic orbits 


In this section we examine how some of the most impor¬ 
tant families of periodic orbits will be influenced by taking 
into account the geometry of a B/P. To do this we study 
two models: the model with the fiducial peanut bulge height 
function, and the model with the sech^ height function with¬ 
out a B/P bulge. We set the pattern speed of both models 
to be such that corotation occurs just outside the bar ra¬ 
dius, within the range 1.4 > Rcn/Rbar > 1, where Rcr 
and Rbar are the corotation and bar radius respectively (e.g. 
Athanassoula|pT92b ). The orbits are calculated in a frame 
of reference co-rot at ing with the pattern speed of the bar. 


In Fig. 6(a)| we show a few typical orbits in the bar 
region for the potentials we are examining, overplotted on 
the image of our simulated galaxy, gtrllG, shown face-on. 
The three most important families of orbits in the bar region 
are shown, i.e. the xi (red lines, extended along the bar 
major axis), X 2 (cyan lines, perpendicular to the bar major 
axis) and 3/1 (green lines, asymmetric with respect to the 
y-axis) families, which are stable along most of their extent. 

In Fig. |6(b)| we plot the characteristic diagram of pe¬ 
riodic orbits, for the two cases with and without a peanut. 
The characteristic diagram gives the value at which the or¬ 
bit intersects the y-axis (yo) as a function of its Jacobian 
energy {Ej, i.e. energy in the rotating frame of reference; 
Binney & Tremaine ( 20Q8| )). The Jacobian energy is in ar¬ 
bitrary units, since, as already mentioned, the mass is also 
in arbitrary units. We see that the characteristic diagram 
of the two models differs significantly. The most noticeable 
effect due to the presence of a B/P is the change in the bifur¬ 
cation loci of the upper and lower branch of the 3/1 family. 
This indicates that taking into account the geometry of a 
B/P in the model changes the location of the 3:1 resonance, 
and therefore the 3/1 family of periodic orbits appears at 
higher energies. Thus orbits of the 3/1 family will differ in 
the two cases, as can be seen in Fig. |7(b)[ where we plot 
two 3/1 orbits in the two models for the same cut along the 
y-axis. The extent of the 3/1 family of orbits is also signif¬ 
icantly increased for the case with a B/P bulge, surpassing 
the extent of the xi family of periodic orbits, which is in fact 
shorter compared to the xi family in the model without a 
B/P bulge. 

The xi family also suffers changes, in the Ej region 
between -1.1 and -0.8. In this area of the diagram, the max¬ 
imum extent of the orbits along the a:-axis reaches the region 
where the effect of the B/P is maximum; therefore, for these 
energies the orbits of the two models differ. In Fig. |7(a) 


show an orbit of the same energy (Ej =-0.9) plotted in the 
two models. When a B/P is present, the maximum extent of 
the orbit along the x-axis is reduced by 12%, while its max¬ 
imum extent along the y-axis is reduced by 46% (measured 
at x=0). 

For the X 2 family, the highest (lowest) value of yo in¬ 
creases (decreases) for the model with a B/P bulge (see 
Fig. |6(b)'] ), i.e. the entire extent of the X 2 family is increased 
by about 43% . As the extent of the X 2 family is related 
to the distance between the two Inner Lindblad Resonances 
(ILRs; [Athanassoula| 1992a ), the distance between these two 
resonances will therefore also increase. This increase is due 
to a weakening of the non-axisymmetric perturbation: when 
the geometry of B/Ps is taken into account in the model, the 
scaleheight of the galaxy is increased, where the B/P is max¬ 
imum, and therefore - since the amount of mass is the same- 
the volume density in the plane of the galaxy is decreased. 
This leads to a decrease in the radial and tangential forces 
in such a way that the non-axisymmetric perturbation is di¬ 
minished, thus changing the distance between the two ILRs. 
This is in accordance with results from both |Contopoulos| 
& Grosbol (1989) and Athanassoula (1992a), the latter of 


which showed that there are a number of model parameters 
which can affect this distance. In particular, in Figures 6 & 


7 of Athanassoula (1992a), we see that the distance between 


the ILRs can increase due to a decrease in the bar mass or 
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(a) Typical orbits in bar region (b) Characteristic Diagram 


Figure 6. Left: Some typical orbits in the bar region for the model with a sech^ height function, over-plotted on the image of the 
simulated galaxy gtrllG: In red the xi bar supporting orbits, in cyan an X 2 orbit perpendicular to the bar and in green the 3/1 orbits 
(asymmetric with respect to the y-axis). Right: Characteristic diagram (intersection of each orbit with the y-axis as a function of the 
Jacobi energy) for the models created from the image of the simulated galaxy gtrllG and the two height functions. The solid black line 
gives the characteristic diagram for the model with the fiducial peanut bulge and the dashed red line the characteristic diagram for a 
model with a fiat sech^ height function. The dotted blue line shows the zero velocity curve (ZVC) for the sech^ model (the ZVC of the 
two models are very similar). 



X (kpc) 

(a) Comparing xi orbits 


X (kpc) 

(b) Comparing 3/1 orbits 


Figure 7. Left: Two xi orbits with the same Jacobian energy (£'j=-0.9) calculated in the two potentials: with (solid black line) and 
without (dashed red line) a B/P bulge. We see that in the B/P model the xi orbit’s length and height are reduced (its extent along the 
aj-axis is reduced by ~12% and along the y-axis by ~4G% -measured respectively at y=0 and a:=0). Right: Two 3/1 periodic orbits for 
yo=1.2 in the two potentials (colours as before). Again the orbits differ significantly. 


the pattern speed, or due to an increase in the central mass 
concentration of the galaxy or the axial ratio of the bar. 

The differences between the orbital families of the two 
models will have effects on their stellar, as well as their 
gaseous dynamics. The extent of the X 2 orbits 
role on the shape of the shock loci in the gas 


plays a crucial 
(Athanassoula 


1992a|b ), and therefore the shape of the gas shocks in the 


two models should differ significantly; conversely, the shape 
and strength of the shocks influence the amount of gas inflow 
towards the centre of the galaxy, and it is therefore likely 
that there will be a measurable difference in the amount of 
gas inflow in models with and without a B/P bulge. This 
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is further supported by the results in Section [4.3| where we 
show that B/P bulges reduce the strength of the bar; we 
plan to address in upcoming work the extent to which the 
gas flows will be affected. 


4.3 B/P effect on bar strength 


We study the effect of the B/P bulge on one of the mea¬ 
sures of bar strength, which involves calculating the non- 
axisymmetric forcings on the disc due to the bar, i.e. the 
bar-induced torque ( Combes Sanders|1981 Buta & Block 


2001). The magnitude of this non-axisymmetric perturba¬ 


tion is given by 


Qrir) = 




(16) 


where Ft is the tangential force FT(r) = (l/r) {d^/dcj)), and 
(Fi^(r)) is the average over azimuth of the radial force 
FR{r)=d^/dr. The forces are calculated directly from the 
image, as described in Section In order to obtain a single 
measure of the bar strength for a galaxy, the quantity Qb, 
which is the maximum of Qt in the bar region, is commonly 
defined as the bar strength. In what follows we investigate 
the effect that a B/P height function will have on both Qb 
and Qt. 

The results of this study are discussed in paragraphs 
(z), (ii) and (m), where we examine models with the flat, 
the fiducial peanut and the fiducial boxy height function re¬ 
spectively. In Table we show the maximum and average 
relative errors of Qt, denoted MAX(Error Qt) and (Error 
Qt) respectively, as well as the relative error of Qb, when 
comparing two models with different height functions. The 
average relative error of Qt over radius is calculated accord¬ 
ing to: 


I "" 

(Error Qt) = — x ^(abs( 


Qrijrj) - Qt2(^0> 

QTi(ri) 


xIOO). (17) 


Plots of these results can also be seen in Eig. where 
it is worth noting that for the flat and peanut height func¬ 
tions, the strength of the bar is Qb~0.55-0.6. According to 


Buta & Block (2001), this represents a strong bar case (be¬ 
tween bar class 5 and 6), which corresponds to approxi¬ 
mately 20% of their sample of SB galaxies. We have already 
shown in Eig. the striking morphological similarity be¬ 
tween IC 4290 and our galaxy, and we note here that IC 4290 
is also classified by |Buta Block (2001) as a class 6 barred 
galaxy, with Q^^O.SG. Therefore the results presented in 
this section, as well as in previous and subsequent sections, 
correspond straightforwardly and quantitatively to strongly 
barred galaxies. However, even weakly barred galaxies will 
have B/P bulges, albeit weaker ones, and therefore the re¬ 
sults will also apply to these galaxies although to a lesser 
extent. We intend to carry out a full statistical study of the 
effects of different strength B/P bulges on the models of their 
host galaxy, together with a full comparison to observations, 
elsewhere. 

(i) A model with a seek and secF height funetion 

In previous work by Laurikainen & Salo (2002) the 


effect of position-independent height functions and height 


Table 1. Errors in bar strength 


Comparison 

(Error Qt) 

MAX(Error Qt) 

Qb 

sech - sech^ 

1% 

1.6% 

1.5% 

peanut - no peanut 

27% 

74% 

4% 

peanut - boxy 

14% 

42% 

16% 


We show the average and maximum of the relative errors of Qt, 
as well as the relative error of Qb, for three different comparisons 
of the setups. 


functions which only vary as a function of radius was 
examined, and these were found not to change Q^ in a 
significant way. Eor realistic height functions such as the 
exponential, sech, or sech^ models, they found that Qb was 
affected by less than 5%, which is consistent with our own 
results. This is confirmed in Eig. and Table where we 
see that for an equivalent scaleheight, the height functions 
of sech and sech^ will produce very similar bar strengths, 
which will tend to become even more similar the thinner 
the disc. 

(ii) A model with the fiducial peanut height function 

We plot Qt for our fiducial B/P height function. In 
the region around the scaleheight maxima (at r=3), Qt is 
significantly flatter than the model without a peanut, due 
to the reduction of the strength of the tangential and radial 
forces. The value of Qb will not be significantly different 
from the case with the flat height function, due to the 
fact that the maxima of the peanut and the maximum 
of Qt are at a relatively large distance from each other. 
The torque induced by the bar in the two cases however is 
significantly different, as can be seen both in the plot and 
in the second row of Table By using the Qb method of 
measuring bar strength, the bars will be judged as having 
the same strength (class 6), and hence the same effect 
on the disc, even though the forces in the plane of the 
galaxy are significantly reduced in the presence of a peanut. 
Therefore, in Section |5.3| we introduce another measure 
of bar strength, which can capture the reduction in bar 
strength when a B/P bulge is present. 

(iii) A model with the fiducial boxy height function 

We also plot Qt for our fiducial boxy height function. 
We see again that where the boxy bulge is maximum, Qt 
is flattened due to the decrease in the strength of the bar 
forces in the plane. We also see that where the top-hat boxy 
function ends, Qt exceeds the values of the sech and sech^ 
curves. This is due to the effect of the boxy bulge on the 
tangential and radial forces, with the former increasing just 
outside the boxy bulge while the latter is decreased in the 
whole disc due to the overall decrease in mass-density in the 
plane of the galaxy. The combination of these two effects 
results in the torque becoming large in the region just out¬ 
side the boxy bulge. As a consequence, Qb is overestimated 
by 16% compared to the fiducial peanut case (third row of 
Table [^. We investigate the boxy height function as an al¬ 
ternative to using the peanut height function - since there is 
one parameter less to model - and conclude that even if one 
is merely interested in Qt, a simple boxy height function is 
not a good approximation to a peanut function. 
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Figure 8. Strength of non-axisymmetric forcings (Qt) as a function of radius, for models with different height functions: sech (solid red 
line), sech^ (dashed green line), the fiducial peanut setup (thick black solid line) and the fiducial boxy setup (dash-dotted blue line). The 
vertical solid black line indicates the radius at which the scaleheight of the fiducial peanut is maximum and the vertical dashed black 
line indicates the end of the bar. 


5 ERRORS DUE TO BOXY/PEANUT 
MODELLING 

In this section we investigate how much error will be induced 
if we include a B/P bulge in the model, but with the main 
peanut parameters differing from that of our fiducial model. 
This type of error is induced due to observational uncertain¬ 
ties, as it is not trivial to observationally obtain the correct 
parameters for the B/P bulge we want to model. This is due 
to the physics of the problem, not the numerical part of the 


calculation (as is the error referred to in Section 2.1 which 


can be made arbitrarily small) and is therefore practically an 
unavoidable source of uncertainty. Nevertheless, as we will 
discuss below, there do exist empirical and theoretical ar¬ 
guments which can constrain the parameter space of a B/P 
bulge. 

The height function we have chosen for our fiducial B/P 
bulge, the sum of two two-dimensional gaussians, has three 
degrees of freedom. Thus inaccuracies in the modelling of the 
B/P bulge can also be introduced in three ways: by estimat¬ 
ing wrongly the height of the gaussians (which corresponds 
to a change in peanut strength), or the distance between the 
maxima of the gaussians (which corresponds to a change in 
the peanut length), or the widths of the gaussians (which 
corresponds to a change in peanut ‘width’, i.e. how peaked 
or thin the peanut is at its maximum). 


5.1 Potential and forces 

In this subsection we investigate how much error is intro¬ 
duced in the potential and forces by incorrectly modelling 
the B/P bulge. 


5.1.1 Peanut strength uneertainties 


The maximum value of the scaleheight of the peanut, also 
called the peanut strength, is a value which is not trivial 
to find observationally. Numerical studies have shown that 
the strength of the peanut correlates with the bar strength 
(Athanassoula 2006). However this relation has a consid¬ 
erable scatter, and can merely give an approximate esti¬ 
mate. Debattista et al.|(2005) showed that for face-on, or 


nearly face-on A-body simulated galaxies, an indicator of 
the presence and strength of a B/P bulge is the fourth or¬ 
der Gauss-Hermite moment of the line-of-sight velocity dis¬ 
persion, /i 4 . However this relation has not been quantified 
in such a way which would allow a direct measurement of 


peanut strength from . Studies of orbital structure (Patsis 
et al.|2002 Skokos et al.|200^ ) have also suggested specific 


families of periodic orbits which are responsible for giving 
the B/P bulge its height, but no direct measurement of the 
strength of the B/P bulge is available from orbital structure 
either. 


In order to measure the error due to peanut strength 
uncertainties, we use a grid of models with varying peanut 
strength and compare them to our fiducial peanut setup. 
These results can be seen in Table [2] below. In this and in 
all subsequent tables, the term ‘Average Error’, indicates 
the average error within the outer isophote of the bar and 
‘Maximum Error’ corresponds to the maximum error found 
in the grid, excluding the central-most point. We see that an 
over- or under-estimation of the error by the same amount 
will produce similar errors in the model. In all the cases 
studied, the error induced by an incorrect peanut strength 
is always less than that induced by not modelling a B/P 
bulge at all. 
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Table 2. Percentage Error due to Peanut Strength Uncertainty 


Peanut Strength 

Average Error 

Maximum Error 


Fx 

Py 


Fx 

Fy 

+50% 

1 . 3 % 

2% 

3% 

3 % 

10 % 

6% 

+25% 

0 . 7 % 

1% 

2% 

1 . 5 % 

5 % 

3% 

-25% 

0 . 7 % 

1% 

2% 

1 . 5 % 

6% 

4% 

-50% 

1 . 4 % 

3% 

4% 

3 % 

14% 

10% 

no peanut 

3 % 

8% 

10% 

7 % 

37% 

28% 


Average and maximum errors of the potential and forces for setups 
with different peanut strength error, within the area enclosed by 
the outer isophote of the bar. The last row gives the error induced 
by not including a B/P bulge in the model at all. 

5.1.2 Peanut width uncertainties 

In Tablewe show the errors for a grid of models with differ¬ 
ent peanut width errors. For all cases considered, the error 
induced due to a miscalculation of the peanut width is less 
than that induced by not modelling a peanut at all, apart 
from the maximum error induced in the potential when the 
peanut width is 50% larger than in the fiducial scenario. This 
is due to a sharp increase in scaleheight in the central re¬ 
gion for large peanut widths (see solid red line in Fig. |lQ((^ , 
which is where the potential is most affected. This error how¬ 
ever is confined only to the potential and to the central most 
grid points, and should not have significant effects on orbital 
calculations in most of the galaxy. 


5.1.3 Peanut length uncertainties 


Of the three parameters - length, strength and width - 
length has the least uncertainty, due to a method proposed 


by Athanassoula et al. (2014). The method determines the 
length of B/Ps for face-on and moderately inclined galaxies, 
which uses the shape of the projected isophotes in the bar 
region. They demonstrated that the barlens and the peanut 
are the same component and therefore that the size of the 
former can be used to estimate the length of the latter. For 
galaxies with larger inclinations the length can be estimated 
from other morphological features in the isophotes created 
by the B/P bulge ( Athanassoula Beaton||2006| Erwin 
Debattista|2013 ), while orbital structure studies confirm the 
aforementioned results and also give clues as to the length 
of the peanut ( Patsis et al.||2002 2003). Due to all this, the 
uncertainties of the length estimates are rather small, cer¬ 
tainly smaller than the corresponding ones for strength and 
width, which is why we use a smaller range of uncertainties 
for the peanut length. 


We carry out comparisons for a grid of models with dif¬ 
ferent peanut length errors and give the results in Table 
As expected, the more we change the length of the peanut 
away from the fiducial value, the larger the errors will be, 
although there is an asymmetry in the error induced with 
respect to over- and under-estimating the length; by under¬ 
estimating the peanut length by a certain amount, we induce 
more error than by overestimating it by the same amount. 
By decreasing the length of the peanut we induce more er¬ 
ror in the central regions of the galaxy, which is where the 
potential is most affected, due to the two gaussians overlap¬ 
ping in the centre and thus increasing the scaleheight (see 
dotted magenta line. Fig. |10(^ ). This can be seen in Table 
1^ for the case of -30% peanut length where, for the poten¬ 
tial, the maximum and average errors are larger than that 
of the +30% peanut length case. 

For all the cases considered, the average error induced 
is smaller or equal to that of not modelling the B/P bulge. 
Given that the length is a fairly well constrained quantity, 
large errors are not expected to be present due to the peanut 
length in the modelling of the B/P bulge. 


5 . 1.4 Combinations of uncertainties 

It is likely that a combination of different kinds of error will 
contribute to the total error budget of a model of the B/P. It 
is not in the scope of this paper to explore the full parameter 
space of the possible error combinations, instead we choose 
a few cases in order to get a feel of the amount of error 
that can be induced. By ‘combination of error’ we refer to a 
combination of all the different sources of error. By ‘+50%’ 
(‘-50%’) we refer to a setup with peanut strength and width 
which are 50% larger (smaller) than the fiducial value, and 
a peanut length which is 30% larger (smaller) than the fidu¬ 
cial. We do not find it necessary to further increase the error 
in peanut length, since it is the best constrained quantity out 
of the three parameters. The ‘+25%’ (‘-25%’) setup corre¬ 
sponds to one with peanut strength, width and length which 
is 25% larger (smaller) than the fiducial value. In Tablewe 
see that all the combinations of uncertainties will introduce 
less error in the model than not including a B/P bulge at 
all. 


5.2 Periodic orbits 


As shown in Section 4.2 the presence of a B/P bulge will 
affect the extent and shape of the different families of peri¬ 
odic orbits which make up the bar. In this section we qual- 


Table 3. Percentage Error due to Peanut Width Uncertainty Table 4. Percentage Error due to Peanut Length Uncertainty 


Peanut Width 

Average Error 

Maximum Error 

Peanut Length 

Average Error 

Maximum Error 


Fx 

Fy 


Fx 

Fy 


Fx 

Fy 


Fx 

Fy 

+50% 

3% 

5% 

6% 

15% 

26% 

27% 

+30% 

1% 

6% 

5% 

4 % 

24% 

13% 

+25% 

1% 

2% 

3% 

4% 

10% 

12% 

+16% 

0.9% 

3% 

3% 

3 % 

12% 

7% 

-25% 

1% 

2% 

3% 

3% 

11% 

6% 

-16% 

1% 

4% 

4% 

5 % 

11% 

10% 

-50% 

2% 

5% 

6% 

5% 

11% 

15% 

-30% 

3% 

8% 

9% 

13% 

25% 

28% 

no peanut 

3% 

8% 

10% 

7% 

37% 

28% 

no peanut 

3% 

8% 

10% 

7% 

37% 

28% 


Average and maximum errors of the potential and forces for setups 
with different peanut width errors within the area enclosed by the 
outer isophote of the bar. The last row gives the error that would 
be present if we do not model a B/P bulge at all. 


Average and maximum errors of the potential and forces for setups 
with different peanut length errors within the area enclosed by the 
outer isophote of the bar. The last row gives the error induced by 
not modelling a B/P bulge at all. 
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Figure 9. (a) Characteristic diagram for models with different B/P setups. See Fig. ^and the text in Section [4.2| for more details on the 
interpretation of the characteristic diagram. The dotted blue line gives the zero velocity curve; the characteristic diagram for the model 
with 0%, 50% and 100% the fiducial peanut strength is given by the dashed red line, the magenta dashed-dotted line and the solid black 
line respectively, (b) The same 3/1 orbit, which cuts the y-axis at y=l, in three models: without a B/P bulge (dashed red line), with 
50% peanut strength (magenta dashed-dotted line) and with 100% peanut strength (black solid line). 




(a) (b) 


Table 5. Percentage Error due to Combination of Uncertainties 


Combination 

Average Error 

Maximum Error 


Pa. 

Fy 

<F 

Pa. 

Fy 

+50% 

2 % 

3% 

5% 

6 % 

23% 

22 % 

+25% 

1 % 

4% 

4% 

6 % 

19% 

18% 

-25% 

1 % 

6 % 

6 % 

5% 

29% 

22 % 

-50% 

2 % 

7% 

8 % 

6 % 

35% 

26% 

no peanut 

3% 

8 % 

10 % 

7% 

37% 

28% 


Average and maximum errors of the the potential and forces for 
setups with different combinations of errors within the area en¬ 
closed by the outer isophote of the bar. The last row gives the 
error induced by not modelling a B/P bulge at all. 


itatively explore the errors introduced in the calculation of 
periodic orbits due to incorrect modelling of a B/P bulge. To 
do this we examine the characteristic diagram of the most 
relevant families of periodic orbits for three models with dif¬ 
ferent peanut strengths (100%, 50% and 0% of the fiducial 
strength), shown in Fig. |9(a)] 

For the model with 50% the fiducial strength the extent 
of the X 2 orbits is reduced by ~19%, while if we do not add a 
B/P at all, the extent of the X 2 family is reduced by ~43% 
(more than double the error for the 50% peanut strength 
case). 

The bifurcation locus of the 3/1 family for the model 
with 50% the peanut strength occurs about halfway between 
the locus of the models with and without a B/P bulge. Ad¬ 
ditionally, the extent of the 3/1 family in the characteristic 
diagram for this model is almost the same as for the model 
with the fiducial peanut, while the extent of the 3/1 fam¬ 
ily without a B/P is significantly shortened. In Fig. |9(b)] we 
can see how the 3/1 orbits are affected by the incorrect mod¬ 
elling of the B/P bulge: for the same cut along the y-axis the 
3/1 orbits are more elongated in the case with the fiducial 
peanut model, while they become less elongated and more 


concave with respect to the bar as the peanut strength is 
reduced. However, as expected, the orbits in the model with 
50% the peanut bulge better match the orbit of the fiducial 
B/P setup than the model without the B/P bulge. 

As already noted in Section [4.2| the xi family of orbits 
is also affected by the presence of the B/P bulge, when the 
maximum extent of the orbits reach the region where the 
effect of the B/P is maximum, i.e. around (x, y)=(zb3kpc, 
Okpc). On the characteristic diagram this occurs in the re¬ 
gion around {Ej, yo)=(0.9, 0.5kpc). However even by un¬ 
derestimating the strength of the B/P by 50%, the xi family 
is quite similar to the xi family of the fiducial B/P case. 

We see that in general, the characteristic diagram of 
the model with 50% the fiducial strength has features which 
are more towards the fiducial peanut model and therefore 
even with such a large error in peanut strength, the charac¬ 
teristic diagram of this model reproduces relatively well the 
characteristic diagram of the fiducial B/P model and cer¬ 
tainly better than the model without a B/P bulge. Similar 
results are found when considering errors in peanut width 
and length, and we therefore conclude that it is preferable 
to include a B/P in the model; the orbital structure of the 
model is significantly affected when a B/P bulge is present, 
and by adding a B/P, even with large errors in its parame¬ 
ters, the periodic orbits reproduce the correct structure more 
closely than when not including a B/P at all. 

5.3 Bar strength 

In this section we examine how both the relative errors of Qt 
and those of its maximum value Qb will be affected by uncer¬ 
tainties in the different parameters of the peanut model. We 
also introduce a new measure of bar strength, ^ which 
takes into consideration the integrated bar-induced torque, 
along the entire range of the bar. We do so because even 
though Qb remains relatively unchanged when adding a B/P 
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Peanut Width 
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Peanut Length 
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Figure 10. Top Row: Left: Values of zq for the fiducial peanut strength (solid thick black line), for 50% larger (solid thin red line), 25% 
larger (dashed green line), 25% less (dotted magenta line), 50% less (dot-dashed cyan line) and 0% peanut strength, i.e. an isothermal 
sheet (thin solid black line). Right: Bar-induced torque Qt as a function of radius for models with the aforementioned height functions 
(respective colours). Middle Row: Left: The values of zq for the fiducial peanut width (solid thick black line), 50% larger (solid thin 
red line), 25% larger (dashed green line), 25% less (dotted magenta line) and 50% less peanut width (dot-dashed cyan line). Right: 
Bar-induced torque Qt, for aforementioned models (respective colours). Bottom Row: Left: The values of zq for models with the fiducial 
peanut length (thick solid black line), 30% longer peanut (dashed red line), 30% shorter peanut (dashed-dotted green line) and a model 
without a peanut (i.e. 0% peanut strength-thin solid black line. Right: Bar-induced torque Qt for the aforementioned models (respective 
colours). In all plots, the vertical lines correspond to the positions of the peanut maxima for each respective height function. 






























Effects of B/P Bulges on Galaxy Models 15 


bulge to the model, Qt over its whole range is significantly 
affected (see for example Fig. |lQ(f) and Table [^, and we 
wish to have a measure of this change with a singlr number. 
The bar strength as defined by is given by 


^bar 

Qip* = — / Qt dr, (18) 

Edisc J 
0 

where Vdisc is the disc scalelength. 

To get a good estimate of the difference of Qt between 
two models over the entire radial range, it is best to carry out 
a point by point comparison, and then consider the radially 
averaged relative Qt error. The relative error of is a 
better proxy for this error than although there are cases 
where the relative error of is small, while the average 
relative error of Qt is much more significant (such as the 
first row of Table or vice-versa (first row of Table . 
Therefore it is possible to have two cases with identical Q^t*, 
but locally different Qt- 


5.3.2 Peanut width uncertainties 


We compare setups with varying peanut widths to our fidu¬ 
cial model. Comparisons for Qt can be seen in Figs. lQ(c)| 
and lQ(d)| and as previously mentioned, the mismatch be¬ 
tween the different models is found when the scaleheights 
of the models are different. The extent of the region of Qt 
which is flattened is reduced when the width of the B/P 
is reduced, as expected. Conversely, when we increase the 
width of the B/P bulge, the area of Qt which is flattened is 
increased. 

Values for the errors in the bar region are given in Table 
We see that errors in peanut width do not induce very 
large errors in the average relative error of Qt, compared 
to the errors induced when not including a B/P bulge. The 
errors induced in Q?^^ are not very large either, although Q^, 
in the case of +50% peanut width, has a relative error larger 
than that of not including a B/P bulge. This again shows 
the importance of carrying out a point by point comparison, 
and a comparison of Q^^, in order to determine the errors 
induced in bar strength due to uncertainties. 


5.3.1 Peanut strength uncertainties 


We see in Figs. 10(a) and 10(b) that Qt increases as we re¬ 
duce the strength of the peanut, and it reaches its maximum 
value when the peanut strength is zero, which corresponds 
to the height function of a flat isothermal sheet. 

The values of the average and maximum relative er¬ 
rors of Qt ((Error Qt) and MAX(Error Qt) respectively), 
as well as the relative error of Qh and can be seen in 
Table (and in all subsequent tables in the following sub¬ 
sections) . We see that when we compare an isothermal sheet 
to the fiducial peanut model the error in Qh is of the order 
of 4%. This is not representative of the large change that the 
average relative error of Qt undergoes (27%). This is due 
to the fact that the maximum of Qt does not change much, 
even though Qt itself is affected by a significant amount 


over its entire range (see Eig. 10(b)). On the other hand, 
the change in which takes into account the whole bar 

region, is more representative of the change in the average 
relative error of Qt (20%). 

In all the cases and for all the measurements of bar 
strength, the error introduced in the model due to uncer¬ 
tainty in peanut strength is not as large as the error intro¬ 
duced when not including a B/P bulge in the model. 


Table 6. Percentage Error of Bar Strength due to Peanut 
Strength Uncertainty 


Peanut 

Strength 

(Error Qt) 

MAX(Error Qt) 

Qb 


+50% 

8 % 

17% 

1.4% 

6 % 

+25% 

4% 

9% 

0.7% 

3% 

-25% 

5% 

11 % 

0.5% 

4% 

-50% 

11 % 

26% 

1 .6% 

7% 

no peanut 

27% 

74% 

4% 

20 % 


The error induced in the bar strength due to different amount of 
error in the peanut strength. We see the effect of these uncertain¬ 
ties on the average and maximum relative error in Qt ((Error 
Qt) and MAX (Error Qt) respectively), as well as on the relative 
errors of Qh and Q^^. 


5.3.3 Peanut length uncertainties 


The results of this study are shown in Eigs. |10(e) a nd|10(f)| 
and in TableSomething worth noting in the Eig. |lQ(f)| is 
that the flattening of Qt occurs at the positions where the 
maxima of the peanut are found (which are indicated by the 
corresponding vertical lines). 

In Table [S] we see the errors induced in the different 
measurements of bar strength due to uncertainties in peanut 
length. Eor the case where the peanut length is 30% larger 
than the fiducial value, Qh has a relative error of 17% com¬ 
pared to the error of 4% in Qh when we do not add a peanut. 
This seems to suggest that it can be counter-productive to 
include a B/P in the model when there exist uncertainties in 
peanut length. However, if we examine the average relative 
error in Qt, we see that the error induced in Qt is in fact 
larger when we do not model a B/P than when we miscalcu¬ 
late its length by +30%. This points once again to the need 
for examining the average errors of Qt and not just Qh, as 
the errors induced in Qh are not representative of the error 
induced in Qt- 

We also see in Eig. 10(f)| that even though Qt of the two 
cases (of fiducial peanut and +30% peanut length) differs 
significantly point by point, the area under the curve for 
the two cases is quite similar. This is reflected in the value 
of the relative error of Qt^^, which only suffers a change of 
around 5% while the average error of Qt suffers a change 


Table 7. Percentage Error of Bar Strength due to Peanut Width 
Uncertainty 


Peanut 

Width 

(Error Qt) 

MAX(Error Qt) 

Qb 


+50% 

6 % 

12 % 

9% 

6 % 

+25% 

4% 

11 % 

5% 

4% 

-25% 

7% 

20 % 

3% 

6 % 

-50% 

16% 

39% 

5% 

13% 

no peanut 

27% 

74% 

4% 

20 % 


As in Table but for errors in peanut width. 
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of 16%. We see therefore that is not always a good 

approximation for the average relative error of Qt- 

The important thing to note is that all the cases con¬ 
sidered induce less error in the average error of Qt than not 
modelling the B/P at all. 


5 . 3.4 Combinations of uncertainties 

As has already been discussed, the most likely scenario is 
that of a combination of different sources of error affecting 
our model. The combinations of errors shown in Table [9] are 
as in Section [5.1.4| We see that the average and maximum 
relative error in Qt for all the combinations is less than that 
of not modelling the B/P at all. The scaleheights and bar 
strength for these models can be seen in Figs. |ll(a) and 
|ll(b)| respectively. 


6 SUMMARY CONCLUSIONS 


In this paper we present the effects of a Boxy/Peanut height 
function on the potential, forces, periodic orbits and bar 
strength of a barred galaxy. We show that such height 
functions significantly affect the results, which consequently 
hints to the effects that a Boxy/Peanut bulge will have on 
its host galaxy. 

We present a method for calculating the potential and 
forces due to the stellar component of a disc galaxy, based 
on a three-dimensional integration of the stellar density dis¬ 
tribution, which can be obtained from images of not too in¬ 
clined galaxies combined with a given height function. The 
method gives robust results for different test cases, as well 
as allowing for any general height function to be used, thus 
allowing for complex density distributions to be modelled. 

We used our code on an image extracted from a iV- 
body+SPH simulation of an isolated galaxy, together with 
two flat, posit ion-independent height functions, and two 
position-dependent height functions. Of the two position- 
dependent height functions, one models a peanut bulge and 
one models a boxy bulge. To create an accurate and phys¬ 
ically motivated fiducial height function for the peanut, 
we shaped and fitted our peanut height function to the 
Boxy/Peanut bulge of the simulated galaxy. 

We found, in accordance with previous results in the 
literature ( Laurikainen Salo||20Q2 ), that for the two flat 
height functions the potential and forces do not vary much, 
provided the setups have equivalent scaleheights. This also 


Table 8. Percentage Error of Bar Strength due to Peanut Length 
Uncertainty 


Peanut 

Length 

(Error Qt) 

MAX(Error Qt) 

Qb 

Qp* 

+30% 

16% 

35% 

17% 

5% 

+ 16% 

9% 

19% 

7% 

3% 

+8% 

5% 

10% 

3% 

1% 

-8% 

4% 

10% 

2% 

1% 

-16% 

9% 

19% 

4% 

3% 

-30% 

16% 

39% 

8% 

8% 

no peanut 

27% 

74% 

4% 

20% 


As in Table but for errors in peanut length. 


Table 9. Percentage Error of Bar Strength due to a Combination 
of Uncertainties 


All Errors 

(Error Qt) 

MAX(Error Qt) 

Qb 

Qip-* 

+50% 

19% 

42% 

47 % 

25% 

+25% 

16% 

35% 

28% 

14% 

-25% 

16% 

50% 

7% 

11% 

-50% 

20% 

63% 

7% 

17% 

no peanut 

27% 

74% 

4% 

20% 


As in Table [G] but for different combinations of uncertainties. 


holds true for the bar strength Qt, which does not change 
much for different flat height functions. 

However, we found that for boxy or peanut height func¬ 
tions the potential and forces vary significantly with respect 
to the case in which a flat height function is used (see Figure 
1^. For the potential, the difference can be up to 7% for an 
extended region within the bar. For the difference can be 
as large as 37%, while for Fy this difference can be as large 
as 28%. We therefore concluded that if a Boxy/Peanut bulge 
is present, one should include it when creating a dynamical 
model of the galaxy. 

To further confirm this result, we examined the effect 
of the Boxy/Peanut bulge on the morphology of the most 
important families of periodic orbits found in barred galax¬ 
ies. We see that by taking into account the B/P geometry 
(i.e. by using our fiducial peanut height function) for a given 
energy, the elongation of the xi orbits - the bar-supporting 
orbits elongated parallel to the bar - is decreased; this ef¬ 
fect is most noticeable for orbits in the region where the 
Boxy/Peanut is maximum (around ±3kpc, see Figure [7(8/ ) 
as expected. By adding a Boxy/Peanut to the model the ex¬ 
tent of the X 2 family - the family of orbits perpendicular to 
the bar - is increased by ~43% (see Figure [^(b)| ), as is the 
distance between the two Inner Lindblad Resonances (ILRs). 
Additionally, the position of the 3:1 resonance is changed; 
the 3/1 family - elongated along the bar and asymmetric 
with respect to the y-axis - appears at larger energies and 
is much more extended in the characteristic diagram (see 
Figure [6(1^ ). All the aforementioned effects will have an im¬ 
pact on the stellar as well as the gaseous kinematics of the 
galaxy. The shape and strength of the shocks in the gas will 
be affected, which in turn affects the amount of gas inflow to 
the central parts of the galaxy. This could have an impact 
on the formation of discy bulges and possibly on the fuel 
reservoir for AGN activity. We plan to investigate in future 
work the extent of the effects of B/P bulges on gas flows in 
galaxies. 

We also studied the effect of the Boxy/Peanut bulge on 
the bar strength, as given by the non-axisymmetric forcings 
due to the bar, Qt- The shape as well as the maximum of 
Qt are significantly affected by taking into account the ge¬ 
ometry of a Boxy/Peanut bulge. We found it useful to define 
a new quantity for measuring bar strength, Q^^, which al¬ 
lows us to extract information about the strength of the bar 
by using its whole extent. The presence of a Boxy/Peanut 
bulge, especially at the points where its scaleheight is maxi¬ 
mum, reduces the bar strength (see Figure]^ which confirms 
that the presence of a Boxy/Peanut bulge reduces the bar 
induced torques. 

Even though taking into account the geometry of 
Boxy/Peanut bulges will affect the model, it is not trivial to 
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Figure 11. Left: Values of zq for different combinations of uncertainties. Right: Strength of non-axisymmetric forcings Qt as a function 
of radius, for the different combinations of uncertainties. In order to not over-clutter the plot, the positions of the peanut maxima are 
given by the vertical arrows, from left to right, for the -50%, -25%, fiducial, +25% and +50% cases. 


obtain their parameters observationally. We therefore exam¬ 
ined how much error would be introduced in the results by 
introducing uncertainties in the Boxy/Peanut parameters. 
Each source of error individually (peanut strength, peanut 
length and peanut width), as well as combinations of the 
different sources of error, induce errors in the results which 
in general are considerably less than those induced by not 
modelling the peanut at all. So, for realistic values of uncer¬ 
tainties in the peanut parameters, the error in including a 
peanut will be less than the error induced by not including 
a peanut in the model. 

The simulated galaxy we chose for this study con¬ 
tains a strong bar, corresponding to bar classes 5 and 6 
from the Buta <V Block (2001) classification. Therefore 
the results of this study can be straightforwardly and 
quantitatively applied to real galaxies with similar bar 
and peanut strength, which account for approximately 
20% of SB galaxies in the local universe. Our results are 
also qualitatively relevant to all barred galaxies in the 
secular evolution phase, although for reduced bar and 
peanut strength the effect of the Boxy/Peanut bulge on 
the model is also reduced. In this work we have pre¬ 
sented an in depth study of the effects of a Boxy/Peanut 
bulge on its galaxy model, focusing on a particular test 
case; we plan to present a quantitative statistical study 
of the effects of these bulges on their host galaxies elsewhere. 
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Effects of B/P Bulges on Galaxy Models 


APPENDIX A: THE ANALYTIC GALAXY 
MODEL 

To model the disc we use a Miyamoto-Nagai density 
( Miyamoto fc Nagai|1975 ) which is defined by the potential- 
density pair, 




f{R,z) = - 


GMd 


^jR^ + ia + y'z^ + b^)^’ 


(Al) 


where the fictitious forces due to rotation are taken into ac¬ 
count. 


Pmn{R, z) — 

PMd \ aR^ + (a + ^ + bP-f (A 2 ) 

47r ) [R2 + (a + V«2+62)2]5/2(^2 + ;,2)3/2 ’ 

where R and z are the cylindrical coordinates, G is 
Newton’s gravitational constant, M the total mass of the 
system and a and b are its characteristic lengths. We set 
the parameters a and 6 to 9 and 1.8 respectively, such that 
we obtain a realistic exponential disc with a scalelength of 
about 3kpc with its mass set to 0.56 times the total mass 
of the system ( Gadotti|| 2011 ). 


The bulge is modelled using a Dehnen sphere (Dehnen 


1993), where we set 7 = 0.5 in order to obtain a cuspy density 


distribution. The potential density pair is given by. 


and 


p(r) = 


vbMb 


Stt y/r(rB + ryG ’ 


(A3) 


. . —2GMb / 

i { r ) - ^^(1 - ( 


r -\- pb 




(A4) 


where is a characteristic radius of the system. The mass 
of the bulge is set to 0.34 the total mass of the model which 
is a typical value for the bulge mass ( Gadotti]| 2011 | . 

The bar is modelled using a Ferrers ellipsoid (Ferrers 
1877), whose density is given by. 


where mr is 


po(l — rryy m ^ 1 

0 m> 1 


(A5) 


2 2 2 
X y z 
= -h -^-. 

o ' ' O • 


/32 


T 


(A 6 ) 


The central density of the bar is given by po, while n 
sets the decrease in bar density as a function of position and 
a, /3 and 7 give the sizes of the three semi-principal axes. 
The mass of the bar is 0.1 times the total mass, which is 
again a typical value for real galaxies ( Gadotti]| 2011 ). The 
bar’s semi-major axis is set to a=5kpc, with an axial ratio 
of a/5=2.5, and we use the inhomogeneous n=l case. 

When integrating orbits in these potentials, we do so 
in the rotating frame of reference of the bar, where the bar 
potential rotates with a pattern speed which is set such that 
corotation occurs just outside the end of the bar within the 
range 1.4 > Rcn/Rbar > 1, where Rcr and Rbar are the 
corotation and bar radius respectively (e.g. Athanassoula 


1992b). The equations of motion in a rotating frame of refer¬ 


ence are taken from Chapter 3 of Binney & Tremaine (2008), 


























